
#include "femlib_config.hh"
#include "point.hh"
#include "fe_evaluators.hh"
#include "fe.hh"

/*
  // NOTE:
  // In order for the following specializations to work correctly, we need a templated version of evaluate_basis and evaluate_basis_all.
  // Unfortunately, virtual template methods are impossible and so until the design is changed, rely on FE_Evaluators to do the job.
  //
#ifdef FADBADPP_GRADIENTS

template <unsigned int DIM, fastfem::ElementType ELEMTYPE, unsigned int ORDER, fastfem::PolynomialFamily PFAMILY, fastfem::FunctionSpace FETYPE, typename T>
void FiniteElement<DIM, ELEMTYPE, ORDER, PFAMILY, FETYPE, T>::evaluate_gradient ( 
                                                                            unsigned int ib, 
                                                                            const Point<SPATIAL_DIMENSION, T>& pt, 
                                                                            Point<SPATIAL_DIMENSION, T>& grad
                                                                           )
{
  typedef typename FEGradientTraits<T>::FADGradient GradientType;

  Point<SPATIAL_DIMENSION, GradientType > pp(pt);

  // Differentiate with respect to each primary variable
  for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
    pp[id].diff(id, SPATIAL_DIMENSION);
  
  // Evaluate function and derivatives
  GradientType f = this->evaluate_basis(ib, pp);
  //GradientType f = Imp::evaluate_basis<GradientType>(ib, pp);
  // f.diff(0, 1) ;
  
  // Value of basis function at the requested point
  //double fval=f.x();
  
  // store the value of derivatives: df/dx into input variable 'grad'
  for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
    grad[id] = f.d(id);

  return ;
}
    
template <unsigned int DIM, fastfem::ElementType ELEMTYPE, unsigned int ORDER, fastfem::PolynomialFamily PFAMILY, fastfem::FunctionSpace FETYPE, typename T>
void FiniteElement<DIM, ELEMTYPE, ORDER, PFAMILY, FETYPE, T>::evaluate_gradient_all ( 
                                                                                const PointList<SPATIAL_DIMENSION, T>& ips,
                                                                                std::vector<Point<SPATIAL_DIMENSION, T > >& grad_vals
                                                                              )
{
  typedef typename FEGradientTraits<T>::FADGradient GradientType;
  
  grad_vals.resize( ips.n_points * FETraits<Imp>::N_Nodes ) ;

  for (unsigned int in = 0, offset=0; in < ips.n_points; ++in, offset+=FETraits<Imp>::N_Nodes )
  {
    Point<SPATIAL_DIMENSION, GradientType > pp(ips.get_point(in));

    // Differentiate with respect to each primary variable
    for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
      pp[id].diff(id, SPATIAL_DIMENSION);

    for (unsigned int ib = 0; ib < FETraits<Imp>::N_Nodes; ++ib)
    {
      // Evaluate function and derivatives
      GradientType f = this->evaluate_basis(ib, pp);
      //GradientType f = Imp::evaluate_basis<GradientType>(ib, pp);
      // f.diff(0, 1) ;
      
      // Value of basis function at the requested point
      //double fval=f.x();
      
      // store the value of derivatives: df/dx into input variable 'grad'
      for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
        grad_vals[offset+ib][id] = f.d(id);
    }
  }

  return ;
}
      
#endif


#ifdef FADBADPP_HESSIANS
        
template <unsigned int DIM, fastfem::ElementType ELEMTYPE, unsigned int ORDER, fastfem::PolynomialFamily PFAMILY, fastfem::FunctionSpace FETYPE, typename T>
void FiniteElement<DIM, ELEMTYPE, ORDER, PFAMILY, FETYPE, T>::evaluate_hessian ( 
                                                                                unsigned int ib, 
                                                                                const Point<SPATIAL_DIMENSION, T>& px, 
                                                                                PointList<SPATIAL_DIMENSION, T>& hess
                                                                              )
{
  femlib_assert(ib < FETraits<Imp>::N_Nodes, "Invalid basis index request");
  
  typedef FEHessianTraits<T>::FADHessian HessianType;
  
  Point<SPATIAL_DIMENSION, HessianType> pp(px);

  // Differentiate with respect to each primary variable
  for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
    pp[id].x().diff(id, SPATIAL_DIMENSION);   // (index id of SPATIAL_DIMENSION)
          
  // Evaluate function and derivatives
  HessianType f = evaluate_basis(ib, pp);
  f.diff(0, 1) ;
  
  // Value of basis function at the requested point
  // double fval=f.x().x();
  // double dfdx=x.d(0).x();    // Value of df/dx
  // double dfdy=y.d(0).x();    // Value of df/dy
  // double dfdxdx=x.d(0).d(0); // Value of df/dx
  // double dfdxdy=x.d(0).d(1); // Value of df/dxdy
  // double dfdydx=y.d(0).d(0); // Value of df/dydx
  // double dfdydy=y.d(0).d(1); // Value of df/dy
  
  // store the value of derivatives: df/dx into input variable 'grad'
  for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
  {
    T* pval = hess.get(id);
    
    for (unsigned int jd=0; jd < SPATIAL_DIMENSION; ++jd)
      pval[jd] = pp[id].d(id).d(jd);
  }
  
  return ;
}

template <unsigned int DIM, fastfem::ElementType ELEMTYPE, unsigned int ORDER, fastfem::PolynomialFamily PFAMILY, fastfem::FunctionSpace FETYPE, typename T>
void FiniteElement<DIM, ELEMTYPE, ORDER, PFAMILY, FETYPE, T>::evaluate_hessian_all ( 
                                                                                const PointList<SPATIAL_DIMENSION, T>& ips, 
                                                                                std::vector<PointList<SPATIAL_DIMENSION, T> >& hess_vals
                                                                              )
{
  typedef FEHessianTraits<T>::FADHessian HessianType;
  
  hess_vals.resize( ips.n_points * FETraits<Imp>::N_Nodes ) ;
  for (unsigned int in = 0, offset=0; in < ips.n_points; ++in, offset+=FETraits<Imp>::N_Nodes )
  {
    Point<SPATIAL_DIMENSION, HessianType > pp(ips.get_point(in));

    // Differentiate with respect to each primary variable
    for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
      pp[id].x().diff(id, SPATIAL_DIMENSION);
 
    
    for (unsigned int ib = 0; ib < FETraits<Imp>::N_Nodes; ++ib)
    {
      // Evaluate function and derivatives
      HessianType f = evaluate_basis(ib, pp);
      f.diff(0, 1) ;

      PointList<SPATIAL_DIMENSION, T>& hess = hess_vals[offset+ib] ;

      // Value of basis function at the requested point
      //double fval=f.x();
      
      // store the value of derivatives: df/dx into input variable 'grad'
      for (unsigned int id=0; id < SPATIAL_DIMENSION; ++id)
      {
        T* pval = hess.get(id);
        
        for (unsigned int jd=0; jd < SPATIAL_DIMENSION; ++jd)
          pval[jd] = pp[id].d(id).d(jd);
      }
    }
  }

  return ;
}
      
#endif

*/

// Some explicit instantiations
// template class FiniteElement<0, fastfem::Vertex, 0, fastfem::Lagrange, fastfem::Pk, double> ;
// template class FiniteElement<1, fastfem::Segment, 1, fastfem::Lagrange, fastfem::Pk, double> ;
// template class FiniteElement<1, fastfem::Segment, 2, fastfem::Lagrange, fastfem::Pk, double> ;

